home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_8.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.0 KB  |  136 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.8 (The Hamming Method).
  9. % Section    9.6,    Predictor-Corrector Method, Page 473
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements Hamming's method
  13.  
  14. % for solving the initial value problem
  15.  
  16. %     y' = f(t,y)    with    y(a) = y0
  17.  
  18. %    Define and store the function f(t,y)    in the M-file  f.m 
  19. % function z = f(t,y)
  20. % z = 30 - 5*y;
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(t,y)');...
  24.            disp('z = 30 - 5*y;');...
  25. diary off;
  26.  
  27. f(0,0); % Remark.  f.m rk4.m hamming.m are usd for Algorithm 9.8
  28. pause      % Press any key to continue.
  29.  
  30. clc;
  31. %    Place the endpoints of [a,b] in  a  and  b.
  32.  
  33. %    Place the initial value y(a) in  ya.
  34.  
  35. %    Place the number of subintervals in  n.
  36.  
  37. a  = 0;
  38. b  = 10;
  39. ya = 1;
  40. n  = 70;
  41.  
  42. pause    % Press any key to continue. 
  43.  
  44. clc;
  45. h  = (b - a)/n;
  46. T = zeros(1,n+1);
  47. Y = zeros(1,n+1);
  48.  
  49. % Generate three starting values using the Runge-Kutta method.
  50. [Ts,Ys] = rk4('f',a,a+3*h,ya,6);
  51. T(1:4) = Ts(1:2:7);
  52. Y(1:4) = Ys(1:2:7);
  53.  
  54. % Proceed with Hamming's method.
  55. [T,Y] = hamming('f',T,Y);
  56. points = [T;Y];
  57.  
  58. pause    % Press any key to see the list of solution points.
  59.  
  60. clc;, clg;
  61. c = 0;
  62. d = 7;
  63. axis([a b c d]);...
  64. plot(T,Y,'-r');...
  65. hold on;...
  66. plot([a b],[0 0],[0 0],[c d]);...
  67. xlabel('t');...
  68. ylabel('y');...
  69. Mx1 = 'Hamming`s solution to y` = f(t,y).';...
  70. title(Mx1);...
  71. grid;...
  72. shg; pause    % Press any key to continue.
  73.  
  74. Mx2 = '     t(k)               y(k)';
  75. Mx3 = 'Hamming`s method is stable';
  76. Mx4 = ['for  n = ',num2str(n),'  and  h = '];
  77. Mx5 = [Mx4,num2str(h),' because'];
  78. Mx6 = '     h < 0.69/|f (t,y)|';
  79. Mx7 = '                y';
  80. clc,echo off,diary output,...
  81. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  82. disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
  83. diary off, echo on
  84.  
  85. pause    % Press any key to continue.
  86.  
  87. clc;
  88. %    Place the endpoints of [a,b] in  a  and  b.
  89.  
  90. %    Place the initial value y(a) in ya.
  91.  
  92. %    Place the number of subintervals in  n.
  93.  
  94. a  = 0;
  95. b  = 10;
  96. ya = 1;
  97. n  = 50;
  98.  
  99. pause    % Press any key to continue. 
  100.  
  101. clc;
  102. h  = (b - a)/n;
  103. T = zeros(1,n+1);
  104. Y = zeros(1,n+1);
  105.  
  106. % Generate three starting values using the Runge-Kutta method.
  107. [Ts,Ys] = rk4('f',a,a+3*h,ya,6);
  108. T(1:4) = Ts(1:2:7);
  109. Y(1:4) = Ys(1:2:7);
  110.  
  111. % Proceed with Hamming's method.
  112. [T,Y] = hamming('f',T,Y);
  113. points = [T;Y];
  114.  
  115. pause    % Press any key to see the list of solution points.
  116.  
  117. clc;
  118. plot(T,Y,'-g');...
  119. xlabel('t');...
  120. ylabel('y');...
  121. title(Mx1);...
  122. grid;...
  123. axis;...
  124. hold off;...
  125. shg; pause    % Press any key to continue.
  126.  
  127. Mx3 = 'Hamming`s method is unstable';
  128. Mx4 = ['for  n = ',num2str(n),'  and  h = '];
  129. Mx5 = [Mx4,num2str(h),' because'];
  130. Mx6 = '     0.69/|f (t,y)| < h';
  131. Mx7 = '            y';
  132. clc,echo off,diary output,...
  133. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  134. disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
  135. diary off, echo on
  136.